%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; clc;
%% manages extraction of data on employment bins
dir = "ThirdChapter_revision";
input = ["/Volumes/Transcend/AWFP/",dir,"/data/CalibrationInput/"];
input = join(input,"");
output= ["/Volumes/Transcend/AWFP/",dir,"/data/CalibrationOutput/"];
output= join(output,"");
cd(input)
ende    = 15;
%%
date = {'23062020'}
cd(input)
d       = csvread('delta_avg_estimation.csv');
N       = exp(csvread('count_avg_estimation.csv'));                 
E       = csvread('employment_ts_estimation.csv');
%% eliminate cities for which delta is smaller than 1:
l = linspace(1,size(d,2),size(d,2));
%%
T = size(E,1);
NC= size(E,2);

DD1      = zeros(1,ende);
DD2      = zeros(1,ende);
select8 = zeros(1,ende);
S7 = zeros(ende,NC);

for x = 1:ende

parametersetting;
delta  = d./(1-alpha);                  % transform the firm size pareto tail into the productivity one
varphi = 1 - 1/(gamma*(1-alpha)+1);
Am     = getAm(N,S,alpha,psi,delta);
Edev_t = E - mean(E,1);
rescale= (varphi./Am).*sqrt(D);         % model backed -> use this one

for t = 1:T-1
   Edevresc(t,:) = Edev_t(t,:)./rescale(t,:);   % this is the X
   Edevresc1(t,:)= Edev_t(t+1,:)./rescale(t,:); % this is the Y
end

for i = 1:NC
   rhom(1,i) = (Edevresc(:,i)'*Edevresc1(:,i))/(Edevresc(:,i)'*Edevresc(:,i));
   varrhom(1,i) = 1/size(E,1)*(Edevresc1(:,i)-rhom(1,i)*Edevresc(:,i))'*(Edevresc1(:,i)-rhom(1,i)*Edevresc(:,i));
end


chi  = psi^(1/(1-alpha));
omega= psi^(-1/(1-alpha));

am = (-varrhom + (rhom-1).*(chi+1)+(1-rhom.^2))./((omega-1)*(chi-omega));
cm = (varrhom - (rhom-1).*(omega+1)-(1-rhom.^2))./((chi-1)*(chi-omega));
bm = 1-am-cm;
deltam  = log(am./cm)./log(psi);

select1 = logical(am<1);
select2 = logical(cm<1);
select3 = logical(am>0);
select4 = logical(cm>0);
select5 = logical(bm<1);
select6 = logical(bm>0);
select9 = logical(deltam>1/(1-alpha));
select7 = select1+select2+select3+select4+select5+select6+select9;

cm   = cm(find(select7==7));
am   = am(find(select7==7));
delta   = delta(find(select7==7)); 
deltam  = deltam(find(select7==7));
select8(1,x) = sum(logical(select7==7));
S7(x,:) = select7;

if size(deltam,2)>1
DD2(1,x) = ((1/select8(1,x))*sum((deltam-delta).^2,2))^.5;
aux = corrcoef([deltam',delta']);
DD1(1,x) = aux(2,1);
else
end

end
plot(select8,DD1,'o')
DD2 = [linspace(1,ende,ende)',DD2',DD1',select8']
DD3 = DD2(find(select8>50),:);
%%
clearvars -except date date1 N E d param input output T NC DD3 sel l; clc;
options = optimset('Display','off');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
loss = [DD3,DD3(:,2) + (NC.*ones(size(DD3,1),1)-DD3(:,4))./(NC.*ones(size(DD3,1),1))];
[~,idx] = sort(loss(:,end)); 
sortedloss = loss(idx,:)  
%x = 5       % select here optimal parameter combination
x = 9     
% change here, a bit arbitrary. but you need to motivate the 
% change in solution to be able to use the new (and correct) data
cd(output)
csvwrite('solution.csv',x) 

parametersetting;
delta  = d./(1-alpha);               % transform the firm size pareto tail into the productivity one
varphi = 1 - 1/(gamma*(1-alpha)+1);
Am     = getAm(N,S,alpha,psi,delta);
Edev_t = E - mean(E,1);
rescale= (varphi./Am).*sqrt(D);    % model backed -> use this one

for t = 1:T-1
   Edevresc(t,:) = Edev_t(t,:)./rescale(t,:);   % this is the X
   Edevresc1(t,:)= Edev_t(t+1,:)./rescale(t,:); % this is the Y
end

for i = 1:NC
   rhom(1,i) = (Edevresc(:,i)'*Edevresc1(:,i))/(Edevresc(:,i)'*Edevresc(:,i));
   varrhom(1,i) = 1/size(E,1)*(Edevresc1(:,i)-rhom(1,i)*Edevresc(:,i))'*(Edevresc1(:,i)-rhom(1,i)*Edevresc(:,i));
end

chi  = psi^(1/(1-alpha));
omega= psi^(-1/(1-alpha));

am = (-varrhom + (rhom-1).*(chi+1)+(1-rhom.^2))./((omega-1)*(chi-omega));
cm = (varrhom - (rhom-1).*(omega+1)-(1-rhom.^2))./((chi-1)*(chi-omega));
bm = 1-am-cm;
deltam  = log(am./cm)./log(psi);

select1 = logical(am<1);
select2 = logical(cm<1);
select3 = logical(am>0);
select4 = logical(cm>0);
select5 = logical(bm<1);
select6 = logical(bm>0);
select9 = logical(deltam>1/(1-alpha));
select7 = select1+select2+select3+select4+select5+select6+select9;

cm   = cm(find(select7==7));
am   = am(find(select7==7));
delta   = delta(find(select7==7)); 
deltam  = deltam(find(select7==7));
rhom   = rhom(find(select7==7)); 
varrhom  = varrhom(find(select7==7));
select8 = sum(logical(select7==7));

S = S(find(select7==7));
N = N(find(select7==7));
l = l(find(select7==7));

Enew = [];
for t = 1:size(E,1)
    aux = E(t,:);
    aux = aux(find(select7==7));
    Enew = [Enew;aux];
end
E = Enew;

if size(deltam,2)>1
DD2 =  (1/select8)*sum((deltam-delta).^2,2);
aux = corrcoef([deltam',delta']);
DD1 = aux(2,1);
else
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NEW %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd(output)
aux = rescale(:,find(select7==7)).*sqrt(varrhom);
csvwrite('instantvolatilitymodel.csv',aux)
csvwrite('instantvolatilitydata.csv',Edev_t(:,find(select7==7)))
csvwrite('selectcities.csv',select7)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OLD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd(input)
toexport = [l',zeros(size(am,2),1),deltam'.*(1-alpha),delta'.*(1-alpha),rhom',zeros(size(am,2),1),am',cm',S']; 
cd(output)
csvwrite('toexport.csv',toexport)

aux = deltam'.*(1-alpha);
toexport = toexport(find(aux>1),:);

c1 = mean(toexport,1);
c2 = var(toexport,1).^.5;
c3 = min(toexport,[],1);
c4 = max(toexport,[],1);

table = [c1(:,2:end);c2(:,2:end);c3(:,2:end);c4(:,2:end)];
csvwrite('tablemoments.csv',table)

table1 = [prctile(S,10),prctile(S,50),prctile(S,90);prctile(am,10),prctile(am,50),prctile(am,90);prctile(bm,10),prctile(bm,50),prctile(bm,90);prctile(cm,10),prctile(cm,50),prctile(cm,90)];
tline = linspace(1,T-1,T-1);
%% 
clear Am Enew Edev_t Edevresc Edevresc1 Emean tline rescale select1 select2 
clear select3 select4 select5 select6 select7 select8 select9 chi omega  
clear D DD1 DD2 toexport bm d

A1 = getAm(N,S,alpha,psi,deltam)./N;
A3 = getAm(N,S,alpha,psi,log(mean(am)./mean(cm))./log(psi))./N;    % variable S
A4 = getAm(N,mean(S),alpha,psi,log(am./mean(cm))./log(psi))./N;    % variable am
A5 = getAm(N,mean(S),alpha,psi,log(mean(am)./cm)./log(psi))./N;    % variable cm

cd(output)
toexport1 = [A1',A3',A4',A5'];
csvwrite('toexport1.csv',toexport1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% find out how well moments for mu match the model implied ones %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:size(A1,2)
    i
    a = am(i);
    c = cm(i);
    b = 1-a-c;
    d = log(a./c)./log(psi);
    p = zeros(S(i),S(i));
    p(1,:) = [a+b,c,zeros(1,S(i)-2)];
    p(2,:) = [a,b,c,zeros(1,S(i)-3)];
        for j=1:S(i)-4
            p(2+j,:) = [zeros(1,j),a,b,c,zeros(1,S(i)-3-j)];
        end
    p(S(i)-1,:) = [zeros(1,S(i)-3),a,b,c];
    p(S(i),:) = [zeros(1,S(i)-2),a,b+c];
    cd(input)
    filename = sprintf('shares%d.csv',i);
    mu = csvread(filename);
    mu = mu';
    pred = p'*mu;
    epsilon = zeros(S(i),size(mu,2));
    varonly = zeros(S(i),size(mu,2));
    VCV = zeros(S(i),S(i),size(mu,2));
    for t = 1:size(mu,2)-1
        v = zeros(S(i),S(i));
    for s = 1:S(i)
            W = diag(p(s,:))-p(s,:)'*p(s,:);
            v = v + mu(s,t).*W;
    end
        epsilon(:,t+1) = mu(:,t+1)-pred(:,t);
        VCV(:,:,t) = v;
        varonly(:,t) = diag(VCV(:,:,t));
    end
    clear VCV v W
    mus = N(1,i)*(1-psi^-d).*psi.^(-d.*linspace(1,S(i),S(i)))/((1-psi^(-d.*S(i))).*psi.^-d);
    musdata = mean(mu,2);
    vnewdata = cov(mu');
    v = zeros(S(i),S(i));
    for s = 1:S(i)
        W = diag(p(s,:))-p(s,:)'*p(s,:);
        v = v + mus(1,s).*W;
    end
    vnew = zeros(S(i),S(i));
    for k = 1:100
        aux = (p')^(k-1)*v*p^(k-1);
        vnew = vnew+aux;
    end
    bin = repmat(linspace(1,S(i),S(i))',size(mu,2),1);
    time = reshape(repmat(linspace(1,size(mu,2),size(mu,2)),size(mu,1),1),size(mu,1)*size(mu,2),1);
    y1 = reshape(mu,size(mu,1)*size(mu,2),1);
    epsilon = reshape(epsilon,size(epsilon,1)*size(epsilon,2),1);
    varonly = reshape(varonly,size(varonly,1)*size(varonly,2),1);
    devmu = repmat(sqrt(mean(musdata-mus).^2)/N(1,i),size(mu,1)*size(mu,2),1);
    devvcv= repmat(sqrt(mean(vnewdata-vnew).^2)/N(1,i),size(mu,1)*size(mu,2),1);
    y2 = reshape(pred,size(pred,1)*size(pred,2),1);
    compact = [time,bin,y1,y2,epsilon,varonly,devmu,devvcv];
    cd(output)
    filename = sprintf('sharesout%d.csv',i);
    csvwrite(filename,compact);
end
%% post-stata
cd(output)
poststata = csvread('r2deltaplot.csv');
plot(poststata(:,1),poststata(:,2),'k.','MarkerSize',10)
lsline
grid on;
ylabel('correlation(\sigma^{data}_{L,m,t},\sigma^{model}_{L,m,t}')
xlabel('\delta^{model}_{m}')
print('-dpng','r2delta')

